
library(rio)
library(ggplot2)
library(texreg)
library(tidyverse)
library(AER)

##set the working directory for the replication archive
working_dir <- "~/Dropbox/dueling project/2nd_paper/replication archive/"

source(paste0(working_dir,"Auxiliary Scripts/RegressionRunner.R"))

##load in the main datasets
newspapers_with <- readRDS(paste0(working_dir,"Data/newspapers_historical.RDS"))
railroads <- readRDS(paste0(working_dir,"Data/rails_waterways.RDS"))
railroads <- 
  railroads %>% 
  arrange(decade) %>% 
  group_by(countyfips) %>% 
  mutate(rail_lag4 = lag(rail_length, 4)) %>% 
  ungroup() %>%
  filter(decade == 1880) %>%
  select(rail_lag1, rail_lag2, rail_lag3, rail_lag4, rail_length, countyfips)

newspapers_with <- inner_join(newspapers_with %>% filter(!is.na(countyfips)), railroads, by = "countyfips")

m0 <- lm(I(all_papers - all_papers1840) ~ log(numPost + 1) + log(total_pop) + log(rail_length + 1) + factor(state), 
         data = newspapers_with)

m1 <- lm(I(all_papers - all_papers1840) ~ log(numPost_lag1 + 1) + log(pop_lag1) + log(rail_lag1 + 1) + factor(state), 
         data = newspapers_with)

m2 <- lm(I(all_papers - all_papers1840) ~ log(numPost_lag2 + 1) + log(pop_lag2) + log(rail_lag2 + 1) + factor(state), 
         data = newspapers_with)

m3 <- lm(I(all_papers - all_papers1840) ~ log(numPost_lag3 + 1) + log(pop_lag3) + log(rail_lag3 + 1) + factor(state), 
         data = newspapers_with)

m4 <- lm(I(all_papers - all_papers1840) ~ log(numPost_lag4 + 1) + log(pop_lag4) + log(rail_lag4 + 1)  + factor(state), 
         data = newspapers_with)

m_base <- lm(I(all_papers - all_papers1840) ~ I(numPost - numPost_lag4) + I(total_pop - pop_lag4) + 
               I(rail_length - rail_lag4) + factor(state), 
             data = newspapers_with)


screenreg(list(m0,m1,m2,m3,m4,m_base), omit.coef = "factor|Median|Hispanic|Gini|Poverty|adults|area",
          custom.model.names = c("Contemporaneous", "First Lags","Second Lags",
                                 "Third Lags","Fourth Lags",
                                 "Change in Num. Newspapers (1884 - 1840)"),
          custom.coef.names =  c("Intercept",
                                 rep(c("Log Post Offices", "Log Population", "Log Rail Length"),5),
                                 "Change in Num. POs (1880-1840)", "Change in Population (1880-1840)", 
                                 "Change in Rail Length (1840-1880)"))

##try making this as a figure
coef_maker <- function(x) coef(x)[2:4]
se_maker <- function(x) summary(x)$coef[2:4,2]

paper_results <- data.frame("Estimate" = c(coef_maker(m0),coef_maker(m1),coef_maker(m2),coef_maker(m3),coef_maker(m4)),
                            "se" = c(se_maker(m0),se_maker(m1),se_maker(m2),se_maker(m3),se_maker(m4)),
                            "Variable" = c(rep(c("Post Offices (logged)", "Population (logged)", "Rail Length (logged)"), 5)), # "Post Office Growth (1840-1880)",
                                           #"Population Growth (1840-1880)"),
                            "Lags" = rev(c(rep(seq(1840,1880,by=10), each = 3))) #, 40, 40))
)

paper_plot <-
  paper_results %>%
  ggplot(. , aes(Lags, Estimate, colour = Variable)) +
  geom_point(size = 5, position=position_dodge(width=3.75)) +
  geom_errorbar(
    aes(ymin = Estimate - 1.65*se, ymax = Estimate + 1.65*se),
    width = 0.01,
    linetype = "solid",
    position=position_dodge(width=3.75)) +
  xlab("Decade of Measurement for Independent Variables") +
  ylab("Effect of Relevant Variable (logged) on\nGrowth in Newspapers (1840-1884)") +
  # coord_flip() +
  theme_bw() +
  theme(legend.position = "bottom", 
        axis.text.x = element_text(size = 25), axis.text.y = element_text(size = 25),
        legend.text = element_text(size = 25),
        axis.title.x = element_text(size = 25),
        axis.title.y = element_text(size = 25),
        strip.text.x = element_text(size = 25)) +
  theme(legend.position="bottom") +
  scale_x_reverse() +
  labs(colour = "") +
  geom_hline(yintercept = 0, linetype = "dashed")

ggsave(paper_plot, 
       file = paste0(working_dir,"Replication Plots/paper_plot.pdf"), 
       width = 12, height = 8)



